function [beta, std_HH] = HHstd_error(y, X, lags)
    % Computes Hansen-Hodrick standard errors
    % y: Dependent variables
    % X: Independent variables
    % lags: Lags to compute HH standard errors
    
    % Estimate Betas and Residuals
    [T,n]   =   size(X);
    beta    =   (inv(X'*X))*X'*y;
    u       =   y-X*beta;
    u= u*ones(1,n);
    err=X.*u;                          % Estimating residuals for each beta

    % HH std errors
    V = [err'*err]/T;                  % Regular weighting matrix
    weight = 0;                        % 0 for even-weighting 1 for Newey-West

    if lags > -1
        for ind_i = (1:lags);
            S = err(1:T - ind_i,:)'*err(1 + ind_i:T,:)/T;
            V = V + (1 - weight*ind_i/(lags + 1))*(S + S');
        end;
    end;
    
    D = inv((X'*X)/T);
    varb = 1/T*D*V*D;
    seb = diag(varb);
    std_HH = sign(seb).*(abs(seb).^0.5);
    
end